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In this paper we demonstrate the direct evidence of soUtons in graphene by means of molecular dynamics 
simulations and mathematical analysis. It shows various solitons emerge in the graphene flakes with two differ- 
ent chiralities by cooling procedures. They are in-plane longitudinal and transverse solitons. Their propagations 
and collisions are studied in details. A soliton solution is derived by making several valid simplifications. We 
hope it shed light on understanding the unusual thermal properties of graphene. 
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Solitons are localized particle-like wavepackets which pre- 
serve their identities such as shape and amplitude during prop- 
agation and after collision between them due to the compen- 
sation of nonlinearity for dispersion. Nowadays, solitons are 
under intense investigation in many systems including Bose- 
Einstein condensates, nonlinear optics, plasmas and anhar- 
monic lattices IJj-^, etc. 

The role solitons played in thermal properties in low di- 
mensional lattices such as thermal rectification and divergence 
of thermal conductivity has been speculated for a long time. 
Thermal conductivity would be direction dependent if solitons 
are involved. It was first theoretically introduced in 1992|6l 
and experimentally observed in asymmetrically mass-loaded 
carbon nanotubes|7| where the rectification coefficient is es- 
timated by considering KdV solitons. Similar thermal recti- 
fication has been numerically studied in asymmetric shaped 
graphene|8| recently. Anomalous thermal conduction that 
thermal conductivity depends on the size of the system has 
been investigated for years ||4| |5l 19^111. Different theoreti- 
cal approaches are applied to understand its mechanism| 12- 
fT4ll . Since solitons appear in those systems with anoma- 
lous transport properties, it is extremely important to assess 
whether such properties are related to solitons. Solitons are 
invoked to explain thermal divergence dates back to Toda in 
1972|[T0l[T5l. Later it shows solitons might be responsible en- 
ergy carriers for thermal divergencel 16] and the divergent ex- 
ponent is dependent upon the scattering rates of solitons and 
soliton-phonon coupled energy diffusion process (Hill. Since 
anomalous thermal conduction was observed in both carbon 
nanotubes and graphene flakes (TTKTSl, thus it is very impor- 
tant to investigate whether solitons exist in graphene. 

Though the thermal properties indicate solitons might exist 
in graphene however there is no direct evidence so far despite 
some theoretical work suggests the possibility of supersonic 
KdV solitons in carbon nanotubes and graphite [19, 20]. Re- 
cently some paper surmises the existence of possible solitons 
in graphene under specific conditions or configurations. Five 
types of solitary waves are numerically found possible along 
the hydrogen-terminated edge of graphene which corresponds 
to the molecule group CH [i21il . When strain is applied spa- 
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Figure 1: (Color online) Chiralities of the rectangle graphene flakes. 
Graphene 6 (0) denotes the graphene flake with armchair edge along 
the X-axis and zigzag edge along the y-axis (orange flake). Since 
rotate Graphene 0(0) clockwise by 30° we get a new rectangle 
graphene flake with armchair edge and zigzag edge swapped (blue 
flake) thus we denote it as Graphene 0(30). 



tially localized electric state might appear near the edge of 
graphene|22|. It is connected to a soliton solution in chiral 
gauge theory which describes the special state of electrons. It 
does not describe the lattice vibration modes. Thus in general 
case whether solitons would exist in graphene is still unan- 
swered. 

In this paper we demonstrate the emergence of solitons in 
graphene in both molecular dynamics simulations and mathe- 
matical analysis. It shows various solitons emerge in graphene 
flakes with two different chiralities by cooling procedures. 
They are in-plane longitudinal and transverse solitons. They 
preserve their identities such as shape and amplitude during 
propagations and after collisions. Their velocities are smaller 
than the relative sound speeds. Phase shifts are brought by an 
averaged acceleration effect in collision between two longitu- 
dinal or two transverse solitons. A longitudinal and a trans- 
verse soliton would simply pass through each other without 
any interactions. We derive a NLS (Nonlinear Schrodinger) 
equation to obtain the analytical soliton solution with several 
simplifications. The validity of the analytical results is dis- 
cussed by comparing with the simulation results. 

We carry out the simulations in two rectangle graphene 
flakes with different chiralities. As shown in Fig. 1, Graphene 
flake with armchair edge along the x-axis and zigzag edge 
along the y-axis is denoted as 6{0) and Graphene flake with 
zigzag edge along the x-axis and armchair edge along the 
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Figure 2: (Color online) Typical solitons in Graphene 0(0) are inves- 
tigated. Above it shows two longitudinal solitons and their velocity 
and displacement distributions. Down it shows their energy and den- 
sity distributions. A pair of small breather- like solitons can also be 
identified. 



y-axis is denoted as 0(30). Graphene 6{0) contains 95952 
atoms and Graphene 0(30) contains 97020 atoms. They are 
485 nm long along the y-axis and 4.7 nm wide along the x- 
axis. Periodic boundary condition is used in the y-axis which 
allows solitons to propagate arbitrary long distance without 
boundary scattering. Periodic or free boundary condition is 
used in the x-axis. Free boundary condition is used in the z- 
axis when concerned. 

We use the second-generation reactive empirical 
bond-order (REBO) potential|23| as implemented in the 
LAMMPS [i24il code to simulate the anharmonic coupling 
between carbon atoms. We take 0.25 fs as the minimum 
timestep. The second-generation reactive empirical bond- 
order (REBO) potential! 23 1 originates from the Brenner's 
potential! 25 1 by including modified fitting data and empirical 
parameters. The general form of both potentials can be 
written as: 



Eij = v/'{nj)^btjV^{nj) 



(1) 



Here V^^ and V^j are repulsive and attractive pairwise po- 
tential of carbon atoms dependent upon their distance, bfj 
is a many-body bond-order term dependent upon their bond 
angles and neighbor atoms. They are widely used to study 
diamond, fullerene, carbon nanotube and graphene, etc. Both 
potentials are capable of describing the short-ranged covalent- 
bonding C-C interactions of the carbon atoms in the same way 
with minor differences. It should be emphasized that, despite 
the effective treatment of hybridization and zero-point energy, 
both potentials include no explicit quantum effects. All con- 
jugation states are derived entirely from the system geometry 



and treat the electronic degrees of freedom in a purely empir- 
ical manner. Thus all the simulation results are obtained and 
discussed in the framework of classic mechanics. 

We also test our simulations by using Tersoff potential !l26l! 
whichi is also a widely used empirical potential for the carbon 
atoms. Similar results would be obtained. It states the main 
results are general for the potentials describing the nonlinear 
interactions of carbon atoms. 

To identify the solitons in graphene a new strategy is used 
in our simulations. Later we shall explain why such strategy 
works. We first thermalize the graphene flakes at a given tem- 
perature (e. g., 1-30 K) for 100 ps with Nose-Hoover thermal 
bath and relax for 20 ps to be equilibrated. Then we set the 
velocities of all atoms to zero every constant step (e. g., 100- 
500 fs) for hundreds of times. After those procedures various 
solitons emerge from graphene with different amplitudes de- 
pendent upon the initial temperature and cooling time. Using 
different initial conditions and time steps, the emerging soli- 
tons might be different. 

We focus upon several different distributions of atoms to 
illustrate the identities of the emerging solitons. Velocity dis- 
tributions of atoms are investigated which include longitudi- 
nal velocity (vy) distributions and transverse velocity (vx) dis- 
tributions. They also stand for the momentum distributions 
when the mass of the carbon atoms is multiplied. Displace- 
ments around the equilibrium positions of the atoms which in- 
clude longitudinal displacement (Y) and transverse displace- 
ment (X) are investigated. Energy distributions are investi- 
gated by considering the differences between the ground state 
energy (AE = E — Eq). Density distributions are investigated 
by considering the reduced density M rather than density p 
for simplicity. The relation between them is: 



p = l.c 



^M- 



M=1.99x 



4V3 



(2) 



Here < r > is the average bond distance of an atom between 
its three neighbors. 

The typical solitons in Graphene 6{0) are shown in Fig. 2- 
4. A pair of longitudinal solitons and a pair of breather-like 
solitons are shown in Fig. 2 and Fig. 3. A pair of longitudi- 
nal solitons and another pair of small longitudinal solitons are 
shown in Fig. 4. 

The velocity and displacement distributions of the two main 
longitudinal solitons are shown in Fig. 2. Their energy and 
density distributions are also shown in Fig. 2. They lack 
energy and aggregate density towards the fluctuation back- 
ground. A pair of small breather-like solitons is also iden- 
tified as shown in Fig. 3. Using different initial conditions 
another pair of small longitudinal solitons can be identified. 
It is shown in Fig. 4 they aggregate energy and lack density 
towards the fluctuation background. 

Typical solitons in Graphene 0(30) are shown in Fig. 5 and 
Fig. 6. A pair of transverse solitons and a pair of longitudinal 
solitons are identified. 

The velocity and displacement distributions of the trans- 
verse and longitudinal solitons are shown in Fig. 5. The 
displacement distributions of the longitudinal solitons are not 
easy to distinguish from the background fluctuation. However 
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Figure 3: (Color online) The breather- like soliton is investigated. 
Above it shows the longitudinal and transverse velocity distributions 
of the breather-like soliton. Down it shows the energy and reduced 
density distributions of the breather- like soliton. 
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Figure 4: (Color online) Using different initial conditions, another 
two small solitons are identified in Graphene d (0) They are also lon- 
gitudinal solitons. 



their shapes resemble the longitudinal solitons in Graphene 
0(0) as well. Their energy and density distributions are shown 
in Fig. 6. The transverse solitons aggregate energy and lack 
density. The longitudinal solitons lack energy towards the 
background. When two transverse solitons collide together 
we can identify them lacking of density then. 

Due to the momentum conservation solitons emerge in pairs 
with opposite velocity localization value. Since the graphene 
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Figure 5: (Color online) Typical solitons in Graphene 0(30) are in- 
vestigated. Above it shows the velocity and displacement distribu- 
tions of the transverse solitons. Down it shows the velocity and dis- 
placement distributions of the longitudinal solitons. 



flakes have no compression or stretching the total displace- 
ments are also zero. The lengths of some solitons are com- 
parable with the mean free path of phononsf27l. It indicates 
a long length might be needed to consider the contribution of 
solitons in thermal conduction. We only observe random fluc- 
tuations in out-of-plane distribution. Such difference might 
affect thermal conduction when flexural modes are enhanced 
or suppressed! 28 1 . 

In order to confirm the wave profiles are solitons we have 
to study their propagation and collision. Solitons shall pre- 
serve their identities in propagation and after collision. It is 
the usual way to identify soliton nature in nonlinear kinetics. 
It is shown in Fig. 7-10 how solitons propagate and collide 
in the graphene flakes. The propagation and collision of the 
longitudinal and breather-like solitons in Graphene d (0) are 
shown in Fig. 7 and Fig. 8. The propagation and collision 
of the transverse and longitudinal solitons in Graphene 0(30) 
are shown in Fig. 9 and Fig. 10. 

It is shown in Fig. 7 the propagation and collision of the 
longitudinal solitons with different amplitudes in Graphene 
0(0). Their propagating velocity is 20.0 km/s and it is smaller 
than the longitudinal sound speed which is 21.3 km/s l|29l . The 
propagation velocity is close to the sound speed. The propa- 
gating velocity seems insensitive to the amplitudes. Their col- 
lisions are elastic without changing their identities. It confirms 
they are solitons. 

It is shown in Fig. 8 the propagation of the breather-like 
solitons in Graphene 0(0). Their propagating velocity is 5.0 
km/s. Their collisions are also elastic. 
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Figure 7: (Color online) Propagation and collision of the two lon- 
gitudinal solitons in Graphene 0(0) with different amplitudes. The 
time interval between two neighboring events is 4.8 ps from (a) to 

(d). 



Figure 6: (Color online) Typical solitons in Graphene 0(30) are in- 
vestigated. Above it shows the energy and density distributions of 
the transverse and longitudinal solitons. Down it shows the energy 
and density distributions when the two longitudinal soliton collide 
together. 



It is shown in Fig. 9 the propagation and collision of the 
transverse solitons in Graphene 0(30). It is shown in Fig. 
10 the propagation and collision of the longitudinal solitons 
in Graphene 0(30). The propagating velocity of transverse 
solitons is 11.4 km/s. The propagating velocity of the lon- 
gitudinal solitons is 20.0 km/s. Both velocities are smaller 
than the relative sound speed which is 21.3 km/s (longitudi- 
nal sound speed) and 13.6 km/s (transverse sound speed) 1291 . 
Their propagating velocities are also weakly dependent upon 
the amplitudes. Their collisions are elastic as well. 

Solitons shall have interactions (decelerations and accel- 
erations) in collisions due to nonlinearity. It is distinctively 
different from linear waves which have no interactions due 
to superposition principle. Solitons would deviate from their 
original trajectories due to collisions. They are not in the po- 
sition which would be anticipated by extending their original 
trajectories. It means solitons are phase-shifted. Phase shift is 
a well known behavior of solitons in collisions 1301 [3T1l. 



It is shown in Fig. 11 the spatial-temporal trajectories of 
the longitudinal solitons in Graphene 0(0). Straight lines 
are observed before or after collision. In collision, solitons 
are first decelerated and then accelerated by a spatial jump. 
The average propagating velocity is 22.3 km/s during colli- 
sion thus solitons exhibit an averaged acceleration effect. The 
dash lines represent the original trajectories of the solitons as- 
suming there was no collision. It clearly illustrates the aver- 
aged acceleration effect brings forth the phase shift of soli- 
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Figure 8: (Color online) Propagation of the breather- like solitons in 
Graphene 0(0). The time interval is 8 ps between two neighboring 
events. 



tons. Spatial jumps are also the important behavior of solitons 
in collision|31|. It is shown in Fig. 11 how the spatial jump 
happens. Since we identify the position of solitons according 
to the extreme value points. So it is difficult to do so in some 
critical moments when there are four extreme value points. It 
behaves like a spatial jump if we neglect those moments. 

It is shown in Fig. 12 the spatial-temporal trajectories 
of two transverse and two longitudinal solitons in Graphene 
0(30). The transverse solitons are first decelerated and then 
accelerated by a spatial jump in collision. The averaged prop- 
agating velocity is 13.9 km/s during collision. It exhibits an 
averaged acceleration effect in collision. The phase shift of 
the transverse solitons is clearly illustrated by the dash lines. 
Similar behaviors are also observed in collision between two 
longitudinal solitons. The averaged propagating velocity is 
22.3 km/s in collision during collision. The phase shift of 
the longitudinal solitons is also clearly illustrated by the dash 
lines. 

It is shown in Fig. 13 the spatial-temporal trajectories of 
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Figure 9: (Color online) Propagation and collision of the transverse 
solitons in Graphene 0(30) with different amplitudes. The time in- 
terval between two neighboring events is 9.6 ps from (a) to (d). 
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Figure 10: (Color online) Propagation and collision of the longitudi- 
nal solitons in Graphene 0(30) with different amplitudes. The time 
interval between two neighboring events is 4.8 ps from (a) to (d). 



a transverse soliton and a longitudinal soliton in Graphene 
0(30). It describes the collision between a transverse soliton 
and a longitudinal soliton. Unlike collision between two trans- 
verse solitons or two longitudinal solitons, they would simply 
pass through each other without any interactions. Thus only 
straight lines are observed in their trajectories. 

Using different initial conditions and cooling time steps 
even more solitons might emerge in graphene with similar ki- 
netic properties. It is shown in Fig. 14 and Fig. 15. 

The emergence of solitons is observed by cooling proce- 
dures. It is quite different from the usual way in which soli- 
tons are generated by imposing a strong impulse upon the 



Figure 1 1 : (Color online) Above it shows the spatial-temporal tra- 
jectories of two longitudinal solitons in Graphene 0(0). Dash lines 
are added as auxiliary lines to illustrate the original trajectories of the 
solitons assuming there was no collision. Down it shows the velocity 
distributions around the spatial jump. The time interval between the 
two neighboring events is 0.0075 ps. ^2 is the critical moment in the 
spatial jump. 



lattice (3TI [32I. Thus a different theoretical approach is re- 
quired to obtain the soliton solution. Here we use the mul- 
tiscale asymptotic method to derive a NLS equation to de- 
scribe the nonlinear process in graphene. Multiscale asymp- 
totic method first developed in 1969 1331 assumes all higher- 
order harmonics are small enough so that the fundamental fre- 
quencies determine the kinetic behaviors of the system. It is 
widely used in different nonlinear systems to obtain soliton 
solutions|34-36|. 

We use the Brenner's potential f25l to describe the interac- 
tions between carbon atoms in the analytical model here. It is 
the 1st simplification in our analytical model. 

Here we take Graphene 0(0) as an example to obtain the 
soliton solution. Following similar steps (191 EOl we expand 
the interatomic potential into Taylor series up to fourth order 
around the equilibrium position in Graphene 0(0). For sim- 
plicity only the longitudinal displacement Y is considered. It 
is the 2nd simplification in our analytical model. 

The Hamiltonian is // = //q + L^ii ^^n where Hn is the per- 
turbation energy in ^th layer. Thus H^ = rrinH^ where m„ is 
the number of atoms in ^th layer and we get 



m 



2(^^+1 



-Ynf 



^(7,+i-y,)3 + ^(F,+i-F,)4(3) 



Here Yn = yn — YnO is the longitudinal displacements of 
atoms around their equilibrium position in ^th layer and b = 
58.93^y, p = 2.796, q = 1.990. 
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Figure 12: (Color online) Above it shows the spatial-temporal tra- 
jectories of two transverse solitons in Graphene 0(30). Down it 
shows the spatial-temporal trajectories of two longitudinal solitons 
in Graphene 0(30). Dash lines are added as auxiliary lines to illus- 
trate the original trajectories of the solitons assuming there was no 
collision. 



Figure 13: (Color online) Above it shows the spatial-temporal tra- 
jectories of a transverse soliton propagate forward and collide with a 
longitudinal soliton in Graphene 0(30). Down it shows the spatial- 
temporal trajectories of a transverse soliton propagate backward and 
collide with a longitudinal soliton in Graphene 0(30). 



Then we can get the kinetic equation to describe the longi- 
tudinal motion as: 



m 



dh 



i [Yn^l ~ 27^ + Yn-\\ — -j2 [\Yn+\ ~ ^n 



bq 



Yn-lY]^Wn^l-YnY -{Yn-Yn-lY] 



{Yn- 

(4) 



Here m = 1.99 x 10"^^kg and k = 1.212A which is the 
equilibrium distance between two layers. 

We then expand the displacement by small excitations up to 
first-order perturbation term 



S^) 



Yn—L,y=l^ ^n \^ ^hniYn) ^ ^^n [T ^ L^ri ^ (pn) 



Wi 



(5) 



Here Un 



.(!)/ 



• Y4L-00 ^n li.'^i in)exp(il(ifn) IS expanded in har- 
monic terms with T = ./-^£^t, t,n = ehIq — Ar/e, 0„ = 

knlo - (OT/e^. 

Here we simplify the expansions again by neglecting the 

high-order harmonics w^i[x^t,n) =0, if |/| > 1. It is the 3rd 
simplification in our analytical model. 

Substituting Yn = e4!o(^'^^) + £u[^^^{T,^ri)exp{i(l)n) + 

£ul^ J*(t, t,n)exp{—i(i>n) into the kinetic equation and equating 
the coefficients of £. After that we derive the following equa- 
tion: 




Figure 14: (Color online) Four solitons in Graphene 0(0). The inter- 
val between the two neighboring events is 12 ps. 




200 300 400 500 

y(nm) 




~§ 0. 



- simulation displacements 

- Fitted displacements 



100 200 300 400 500 

y(mn) 




- simulation velocities 

- fitted velocities 




100 200 300 400 500 

y(mn) 



Figure 15: (Color online) Eight solitons in Graphene 0(30). The 
interval between the two neighboring events is 10 ps. 
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(6) 



Here C is an integration constant. By substituting variables 
it is reduced to a standard cubic NLS equation: 



'l^-fS + lv^lV=o 



'dri 2 35 



(7) 



Here ij/ = exp{iSp^co f CdT)ul j, G 



(p2_ 3^)0)2+4^2' ^ 

cot/g, 5 = ^^n- Here we have G = 433^^313 > which 
means it is a normal NLS equation. 

A soliton solution is possible in the NLS equation above: 



Yn =A[cosX — 2p^]tanh(^ — AtancpsinX + A' 



(8) 



Here A' is an integration constant and A, a, 9 are free pa- 
rameters, vi is the propagation velocity of the soliton. 



knk + 2a - [lacos^ + (2 - a- 
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yn =Yn^YnO = Yn ^ hIq 



Figure 16: (Color online) Comparing the analytical soliton solution 
with the soliton obtained in simulations. Above it shows the displace- 
ment distribution agrees with the simulation results. Down it shows 
the velocity distribution also agrees with the simulation results. 



The complete form of the soliton solution is complicated. It 
can be further simplified by considering the kinetic properties 
of solitons in the molecular simulations. From the previous 
study, we know the velocities of the solitons are insensitive to 
the amplitudes. Thus the amplitude dependent parameter in 
vi might be neglected. It means tan(p = here. It is the 4th 
simplification in our analytical model. 

The propagating velocities of solitons are close to the rela- 
tive sound speed. It also can be used as a valid simplification. 

The relative sound speed is V2 = a/^. It means vi = V2. It is 
the 5th simplification in our analytical model. 
From the 5th simplification we can obtain 



A 



AtarKp ^^.^^klp- 



['^-^?-(« + ff>-f]V^^V^ (9) 

Taking the 4th simplification into account, we have 

1 (10) 

So we may obtain the following relation that: 



A 



A 



cos 



klo 



l,^m^^O, a<Cl 



(11) 



Taking tancp = 0,cos'^^l, sin^ ^ and a <Cl into Eq. 
8, we finally can obtain the soliton solution in the continuum 
approach as: 



Y = Aotanh[c{y - tvi )] + A' 



(12) 



Here Aq and c are determined by the anharmonic parame- 
ters p and q in Eq. 3. 



Ao = A(Vl-4a2 - 2p^/a) = A{Vl-4a^ - 

(13) 



_4p^ 



(p2- 3^)0)2+4^2) 



c=A/v^ = A./(/72-|^)(02 + V 



(14) 



The relative longitudinal velocity distribution of the soliton 
solution is: 

vy=^ = -Aocvisech^[c{y-tvi)] (15) 

In a given time to the relative soliton solution should be: 



Y = Aotanh[c{y — yo)] +A' 
vy = -Aocvisech^[c{y -yo)] 



(16) 
(17) 



To investigate the validity of our analytical results, we first 
try to fit the soliton obtained in the simulation with the soliton 
solution we derived. Next we try to use the fitted soliton so- 
lutions as the initial conditions to investigate the propagation 
and collision of the fitted solitons. 

We first fit the displacement distribution to determine Aq, c 
and A\ When the parameters are obtained we use them to fit 
the velocity distribution. The fitted displacement and velocity 
distributions (vi =20.0 km/s) are: 



Y = -0.044ra^/z[0.04(3;- 100)] +0.044 
vy = 3.5sech^[0.04{y-l00)] 



(18) 
(19) 



It is shown in Fig. 16 the results fit the left side soliton 
well. The right side soliton can also be well fitted in the same 
way. Next we use the fitted solitons as the initial conditions 
to investigate their propagation and collision properties. It is 
shown in Fig. 17 that the propagation and collision proper- 
ties of the fitted solitons are almost the same as the solitons 
obtained in the previous molecular dynamics simulations. It 
states the derived soliton solution can describe the solitons ob- 
tained in the molecular dynamics simulations well . 

It is well known there are extreme difficulties to obtain a 
soliton solution without making proper simplifications. Using 
different methods and simplifications, different soliton solu- 
tions might be obtained. To derive the soliton solution we 
use five simplifications. Thus it is important to understand the 
conditions our simplifications correspond to and the limits of 
these simplifications. 

The 1st simplification is to use the Brenner's potential in- 
stead of the second-generation reactive empirical bond-order 
(REBO) potential. As we point out both potentials describe 
the short-ranged covalent-bonding C-C interactions of the car- 
bon atoms in the same way. Our simulation results are general 
for the potentials describing the C-C interactions. Thus this 
simplification brings no significant differences. 
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Figure 17: (Color online) ^i is the solitons obtained in the previ- 
ous molecular dynamics simulations by cooling procedures. 52 is the 
solitons derived in the analytical results. The propagation and colli- 
sion properties resemble each other well. 




Figure 18: (Color online) Long time behavior of the longitudinal 
solitons in Graphene 0(0). The time interval between two neighbor- 
ing events is 1.2 ns from (a) to (d). 



The 2nd simplification is to consider the displacements 
only in one direction. It is due to the fact that only trans- 
verse/longitudinal solitons would interact with each other 
while a transverse soliton would not interact with a longitu- 
dinal soliton. It simplifies the graphene flakes as a quasi-one 
dimensional chain. However in the molecular dynamics sim- 
ulations, displacements of all directions are considered. Thus 
our analytical model is insufficient to describe the breather- 
like solitons (Fig. 3) obtained in the simulations. In order 
to uncover the origin of the breather-like soliton, a coupled 
NLS equation considering both transverse and longitudinal 
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Figure 19: (Color online) Long time behavior of the transverse soli- 
tons in Graphene 0(30). The time interval between two neighboring 
events is 1.2 ns from (a) to (d). 
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displacements might be needed. 

The 3rd simplification is to neglect the high-order harmon- 
ics in graphene. It is due to the fact that the high-frequency 
modes are suppressed in the molecular dynamics simulations 
by the cooling procedures. Here it also explains why such 
strategy works. However high-order harmonics may take part 
in the solitons obtained in the simulations as well. Thus the 
single soliton solution might be just an approximation for a 
multiple soliton solution. If it is so after a long time the soli- 
tons may split into several secondary solitons if their velocities 
are amplitude dependent. It brings us to the 4th simplification. 

The 4th simplification is to neglect the velocity amplitude 
dependent parameter. It is a valid simplification according to 
Fig. 7, Fig. 9 and Fig. 10. But should the velocities of the 
solitons be absolutely amplitude independent or weakly am- 
plitude dependent still remain unanswered. To answer that 
question we need to investigate the long time behaviors of the 
solitons. Thus we keep the solitons propagate as long as about 
4 ns. The total propagation distance of the longitudinal soli- 
tons is more than 0.1 mm. The total propagation distance of 
the transverse solitons is more than 0.05 mm. It is shown in 
Fig. 18 the longtime behavior of the longitudinal solitons in 
Graphene 0(0). It is shown in Fig. 19 the longtime behavior 
of the transverse solitons in Graphene 0(30). They gradually 
split into several secondary solitons arranged in a train with 
increasing or decreasing amplitudes. If the velocities are pro- 
portional to the amplitudes they would appear with increasing 
amplitudes arrayed in a straight line. Meanwhile if the veloc- 
ities are inversely proportional to the amplitudes they would 
appear with decreasing amplitudes arrayed in a straight line. 
So it is shown in Fig. 20 the velocities of the longitudinal soli- 
tons are proportional and the transverse solitons are inversely 
proportional to the amplitudes. 

The 5th simplification is to assume the soliton velocity is 
close to the relative sound speed. The propagating velocity 



Figure 20: (Color online) Above it shows the longitudinal solitons 
are arranged in a train with increasing amplitudes. Their amplitudes 
appear in a straight line. Down it shows the transverse solitons are 
arranged in a train with decreasing amplitudes. Their amplitudes also 
appear in a straight line. 



of the longitudinal solitons is 20.0 km/s and the propagating 
velocity of the transverse solitons is 1 1.4 km/s. Both velocities 
are close to the relative sound speed which is 21.3 km/s and 
13.6 km/s. 

In summary, various solitons would emerge in the graphene 
flakes by cooling procedures. They are subsonic transverse 
and longitudinal solitons. They exhibit an averaged acceler- 
ation effect in collision which means they are phase-shifted. 
A soliton solution is derived by taking several simplifications 
which correspond to the simulation conditions. The validity 
of the soliton solution is confirmed by comparing with sim- 
ulation results. When periodic boundary condition is used in 
X-axis, the graphene flakes resemble simplified carbon nan- 
otubes. Similar solitons might be also important in carbon 
nano tubes. We hope our work shed light on understanding 
the unusual thermal properties of graphene and other similar 
materials and also benefit the development of graphene based 
thermal device. 
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